////////////////////////////////////////////////////////////////////////////////
//
// labFisInterfElab.cc - Fit data from SNS Lab Fis Interference experiment
// 
// 1st release: A. Lusiani <a.lusiani@sns.it> April 2003
//
// this is a standalone program that does the same as labFisInterfFit.C
//

#include <stdio.h>

#include "TStyle.h"
#include "TApplication.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TH1.h"
#include "TProfile.h"

int main(int argc, char *argv[])
{
  //
  // - open the Root file with data saved from interference experiment
  // - get the microphone voltage vs. time experiment
  // - create an empty monodimensional histogram
  // - fill the new histogram with to voltages vs time from the first one
  // - set the errors in the new histogram
  // - fit a range of the histogram with a sine function
  //
  
  //--- init an object to get Root GUI on a canvas
  TApplication theApp("App", &argc, argv);

  if (!fopen("labFisInterfSave.root","r")) {
    fprintf(stderr, "Cannot find file labFisInterfSave.root!\n");
    return 1;
  }
  TFile* f = new TFile("labFisInterfSave.root");
  
  TProfile *hsMicroTimeProf = (TProfile *)f->Get("hsMicroTimeProf");
  
  UInt_t bins = hsMicroTimeProf->GetNbinsX();
  Double_t xlow = hsMicroTimeProf->GetXaxis()->GetXmin();
  Double_t xupp = hsMicroTimeProf->GetXaxis()->GetXmax();
  
  TH1F *hs = new TH1F("hs", "Microphone voltage vs time", bins, xlow, xupp);
  
  //
  // errors
  //
  // here we are fitting a series of ADC measurements
  //
  // if there is no noise the error of one ADC measurement can be set
  // as 1/sqrt(12) ADC counts (the sigma of a uniform distribution between
  // x and x+1
  //
  // if there is noise, one can repeat the the ADC measurement of a constant
  // voltage and measure the sigma of the ADC count distribution: this
  // is a good estimate of the error to be assigned
  //
  // Here we assume that the error is 4 ADC counts
  // For the interference experiment the ADC range is -5V .. + 5V
  // and the ADC counts are from 0 to 0xffff, according to this the
  // ADC count error is translated to voltage error
  //
  
  Double_t error = 4. * 10. / 0xffff;
  
  //--- build histogram to fit, with data points and errors
  for(UInt_t i=1; i< bins+1; i++) {
    //--- set bin i content (BEWARE: does not correspond to filling)
    hs->SetBinContent(i, hsMicroTimeProf->GetBinContent(i));
    hs->SetBinError(i, error);
  }
  
  //--- fit between 0.00425s and 0.006s
  Double_t fitLow = 0.00425;
  Double_t fitUpp = 0.00600;
  
  // 
  // fit function = offset + A * sin( 2*pi*freq * x + phi )
  // parameter offset accounts for constant voltage offset
  //
  
  TF1 *fun = new TF1 ("fun", "[0]+[1]*sin(atan(1)*8*[2]*x+[3])", fitLow, fitUpp);
  fun->SetParNames( "offset","A", "freq", "#phi");
  //--- initialize parameters to sensible values
  fun->SetParameter( 0, 0);
  fun->SetParameter( 1, 0.05);
  fun->SetParameter( 2, 40000);
  fun->SetParameter( 3, 0);
  
  //
  // set proper limits of parameter A,
  // otherwise the fit does not converge well
  //
  // one could scan the data beforehand to get proper limits
  //
  
  //
  // since the fit has trouble converging, a preliminary fit is done
  // on just few points, in the range 0.0425s - 0.0435s
  //
  
  fun->SetRange(0.00425, 0.00435);
  
  hs->Fit("fun", "r");
  
  //
  // the above preliminary fit converges well and finds good starting point
  // for the parameters. The fit is repeated on the extended range.
  //
  
  fun->SetRange(0.00425, 0.00600);
  
  hs->Fit("fun", "r");

  gStyle->SetOptFit(111);
  hs->Draw();
  
  printf("\n");
  printf("Fit results for f(x) = offset + A * sin( 2*pi*freq * x + phi\n");
  printf("chi2            %6.4f for %d degrees of freedom\n",
       fun->GetChisquare(), fun->GetNDF());
  printf("constant offset %10.4f +- %6.4f mV\n",
	 1000*fun->GetParameter(0), 1000*fun->GetParError(0));
  printf("amplitude       %10.4f +- %6.4f mV\n",
	 1000*fun->GetParameter(1), 1000*fun->GetParError(1));
  printf("frequency       %10.4f +- %6.4f Hz\n",
	 fun->GetParameter(2), fun->GetParError(2));
  printf("phase           %10.4f +- %6.4f\n",
	 fun->GetParameter(3), fun->GetParError(3));
  
  //--- start root GUI: "Exit ROOT" will return to the next statement
  theApp.Run(kTRUE);
}
